# Dickstein, Ho, and Mark (2023)
# This script runs the necesary counterfactuals
# Specifically runs the counterfactuals for counterfactuals involving changing values of psi, 
# so alpha-psi has to change accordingly to an upper boundary on psi

if(!exists("run_which_ctfls")){
  run_which_ctfls <- c(0, 1, 2, 3, 4, 5)
}
if(!exists("create_which_columns")){
  create_which_columns <- c(0, 1, 2, 3, 4, 5)
}

# Base Scenario: No moversif
if(0 %in% run_which_ctfls){
  OnlyIndividualMarketCounterfactual <- RunCounterfactual_psi(
    ExplDat = ExplodedData_Ind_small,
    PlanDat = PlanData,
    behaviorShare = beshr,
    maxR = 150, 
    subsidytypeInd = "Ratio", 
    subsidytypeSG = "Ratio",
    fixed_markup = 0,
    psi_upperbound = psi_upperbound)
}

if(0 %in% create_which_columns){
  resbase <- generateOutcomeColumns_psi(
    CounterFactualResults = OnlyIndividualMarketCounterfactual,
    SGpremsubs = T, 
    subsidytypeInd = "Ratio", 
    subsidytypeSG = "Ratio")
}

#1. Scenario One: they move in to individual market and get individual market subsidies:
if(1 %in% run_which_ctfls){
  CombinedMarketCounterfactual1 <- RunCounterfactual_psi(
    ExplDat = rbind(ExplodedData_Ind_small, ExplodedData_SG_small),
    PlanDat = PlanData,
    behaviorShare = beshr,
    maxR = 150, 
    subsidytypeInd = "Ratio", 
    subsidytypeSG = "Ratio",
    fixed_markup = mrkup,
    psi_upperbound = psi_upperbound)
}

if(1 %in% create_which_columns){
  res1 <- generateOutcomeColumns_psi(
    CounterFactualResults = CombinedMarketCounterfactual1,
    SGpremsubs = T, 
    subsidytypeInd = "Ratio", 
    subsidytypeSG = "Ratio")
}

if(1 %in% create_which_columns & 0 %in% create_which_columns){
  resbase[[2]] <- generateSGOutcomeColumnsPre_psi(
    CounterFactualResults = CombinedMarketCounterfactual1,
    SGpremsubs = T, 
    subsidytypeInd = "Ratio", 
    subsidytypeSG = "Ratio")
}

#2. Scenario Three: ICHRA extension - proportional subsidies
ExplodedData_SG_3 <- ExplodedData_SG
ExplodedData_SG_3[, subsidyshare := 1 - ((1 - (best_guess_sumrate / 100)) * (1 - nu)), ]
ExplodedData_SG_3$silver_sub <- 0
ExplodedData_SG_3$silver_costshare <- 0
ExplodedData_SG_3$x_av <- ifelse(ExplodedData_SG_3$x_av %in% c(0.73, 0.87, 0.94), 0.7, ExplodedData_SG_3$x_av)
ExplodedData_SG_3$fixedeffect_beta0 <- as.matrix(ExplodedData_SG_3[, names(SGNamedBetaList$beta0), with = F]) %*% SGNamedBetaList$beta0
ExplodedData_SG_small_3 <- ExplodedData_SG_3[, .(
  orig_subs_id,
  year,
  hhid = -hhid,
  plan_id,
  best_guess_AV_old,
  best_guess_netprem_old,
  best_guess_sumrate,
  grossprem_old,
  sum_concurrent_risk,
  ratingfactor = RatioSum,
  subs = ifelse(plan_id == "Outside_Option", 0, subsidyshare),
  x_av,
  alpha,
  omega,
  psi,
  fixedeffect_beta0,
  kappa = ifelse(plan_id == "Outside_Option", 1, plan_AV /(x_av * 100))), ]

if(2 %in% run_which_ctfls){
  CombinedMarketCounterfactual3 <- RunCounterfactual_psi(
    ExplDat = rbind(ExplodedData_Ind_small, ExplodedData_SG_small_3),
    PlanDat = PlanData,
    behaviorShare = beshr,
    maxR = 150,
    subsidytypeInd = "Ratio", 
    subsidytypeSG = "Ratio",
    fixed_markup = mrkup,
    psi_upperbound = psi_upperbound)
}

if(2 %in% create_which_columns){
  res2 <- generateOutcomeColumns_psi(
    CounterFactualResults = CombinedMarketCounterfactual3,
    SGpremsubs = F, 
    subsidytypeInd = "Ratio", 
    subsidytypeSG = "Ratio")
}

rm(ExplodedData_SG_3)
rm(ExplodedData_SG_small_3)


#3. Scenario Six: ICHRA extension - lump sum subsidy
ExplodedData_SG_6 <- ExplodedData_SG
ExplodedData_SG_6$silver_sub <- 0
ExplodedData_SG_6$silver_costshare <- 0
ExplodedData_SG_6$x_av <- ifelse(ExplodedData_SG_6$x_av %in% c(0.73, 0.87, 0.94), 0.7, ExplodedData_SG_6$x_av)
ExplodedData_SG_6$fixedeffect_beta0 <- as.matrix(ExplodedData_SG_6[, names(SGNamedBetaList$beta0), with = F]) %*% SGNamedBetaList$beta0

# For those without metal tiers before, estimate amount predicted by characteristics: 
SubDatCollapse <- ExplodedData_SG[, .(
  alpha = mode1(alpha), 
  ratingfactor = mode1(RatioSum), 
  portland = mode1(portland), 
  subsidy = mode1(
    ifelse(is.na(grossprem_old) | is.na(best_guess_sumrate), 
    -999, 
    nu * (1 - (best_guess_sumrate / 100)) * grossprem_old/ 100))), 
  by = c("orig_subs_id", "year")]
lm1 <- lm(data = SubDatCollapse[subsidy > 0], 
          subsidy ~ alpha + portland + ratingfactor + I(ratingfactor ^ 2) + I(ratingfactor ^3))
SubDatCollapse$predictedsubsidy <- predict(lm1, SubDatCollapse)

ExplodedData_SG_6 <- merge(ExplodedData_SG_6, 
                           SubDatCollapse[, .(predictedsubsidy, orig_subs_id, year), ], 
                           by = c("orig_subs_id", "year"), all.x = T)

# Cleaning 
ExplodedData_SG_small_6 <- ExplodedData_SG_6[, .(
  orig_subs_id,
  year,
  hhid = -hhid,
  plan_id,
  best_guess_AV_old,
  best_guess_netprem_old,
  best_guess_sumrate,
  grossprem_old,
  sum_concurrent_risk,
  ratingfactor = RatioSum,
  subs = ifelse(plan_id == "Outside_Option", 0, ifelse(is.na(grossprem_old) | is.na(best_guess_sumrate), predictedsubsidy, nu * (1 - (best_guess_sumrate / 100)) * grossprem_old / 100)),
  x_av,
  alpha,
  omega,
  psi,
  fixedeffect_beta0,
  kappa = ifelse(plan_id == "Outside_Option", 1, plan_AV /(x_av * 100))), ]

if(3 %in% run_which_ctfls){
  CombinedMarketCounterfactual6<- RunCounterfactual_psi(
    ExplDat = rbind(ExplodedData_Ind_small, ExplodedData_SG_small_6),
    PlanDat = PlanData,
    behaviorShare = beshr,
    maxR = 150, 
    subsidytypeInd = "Ratio", 
    subsidytypeSG = "Fixed",
    fixed_markup = mrkup,
    psi_upperbound = psi_upperbound)
}

if(3 %in% create_which_columns){
  res3 <- generateOutcomeColumns_psi(
    CounterFactualResults = CombinedMarketCounterfactual6,
    SGpremsubs = F, 
    subsidytypeInd = "Ratio", 
    subsidytypeSG = "Fixed")
}

# I dont remove scenario 6 datasets, as I use them later. 

#4. Scenario Four: proportional subsidies, forced into insurance
ExplodedData_SG_7 <- ExplodedData_SG
ExplodedData_SG_7[, subsidyshare := 1 - ((1 - (best_guess_sumrate / 100)) * (1 - nu)), ]
ExplodedData_SG_7$silver_sub <- 0
ExplodedData_SG_7$silver_costshare <- 0
ExplodedData_SG_7$x_av <- ifelse(ExplodedData_SG_7$x_av %in% c(0.73, 0.87, 0.94), 0.7, ExplodedData_SG_7$x_av)
ExplodedData_SG_7$fixedeffect_beta0 <- as.matrix(ExplodedData_SG_7[, names(SGNamedBetaList$beta0), with = F]) %*% SGNamedBetaList$beta0
ExplodedData_SG_small_7 <- ExplodedData_SG_7[, .(
  orig_subs_id,
  year,
  hhid = -hhid,
  plan_id,
  best_guess_AV_old,
  best_guess_netprem_old,
  best_guess_sumrate,
  grossprem_old,
  sum_concurrent_risk,
  ratingfactor = RatioSum,
  subs = ifelse(plan_id == "Outside_Option", 0, subsidyshare),
  x_av,
  alpha,
  omega,
  psi,
  fixedeffect_beta0,
  kappa = ifelse(plan_id == "Outside_Option", 1, plan_AV /(x_av * 100))), ]
ExplodedData_SG_small_7 <- ExplodedData_SG_small_7[plan_id != "Outside_Option"]

if(4 %in% run_which_ctfls){
  CombinedMarketCounterfactual4 <- RunCounterfactual_psi(
    ExplDat = rbind(ExplodedData_Ind_small, ExplodedData_SG_small_7),
    PlanDat = PlanData,
    behaviorShare = beshr,
    maxR = 150,  
    subsidytypeInd = "Ratio", 
    subsidytypeSG = "Ratio",
    fixed_markup = mrkup,
    psi_upperbound = psi_upperbound)
}

if(4 %in% create_which_columns){
  res4 <- generateOutcomeColumns_psi(
    CounterFactualResults = CombinedMarketCounterfactual4,
    SGpremsubs = F, 
    subsidytypeInd = "Ratio", 
    subsidytypeSG = "Ratio")
}

# 5. Scenario Seven: forced into insurance, with lump sum
ExplodedData_SG_small_7 <- ExplodedData_SG_small_6[plan_id != "Outside_Option"]

if(5 %in% run_which_ctfls){
  CombinedMarketCounterfactual7 <- RunCounterfactual_psi(
    ExplDat = rbind(ExplodedData_Ind_small, ExplodedData_SG_small_7),
    PlanDat = PlanData,
    behaviorShare = beshr,
    maxR = 150,  
    subsidytypeInd = "Ratio", 
    subsidytypeSG = "Fixed",
    fixed_markup = mrkup,
    psi_upperbound = psi_upperbound)
}

if(5 %in% create_which_columns){
  res5 <- generateOutcomeColumns_psi(
    CounterFactualResults = CombinedMarketCounterfactual7,
    SGpremsubs = F, 
    subsidytypeInd = "Ratio", 
    subsidytypeSG = "Fixed")
}
